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We study a variant of the univariate approximate GCD problem, where the coeffi- 
cients of one polynomial /(a;)are known exactly, whereas the coefficients of the second 
polynomial g{x)ma.y be perturbed. Our approach relies on the properties of the matrix 
which describes the operator of multiplication by gin the quotient ring C[x]/(f). In 
particular, the structure of the null space of the multiplication matrix contains all the 
essential information about GCD(/, g). Moreover, the multiplication matrix exhibits 
a displacement structure that allows us to design a fast algorithm for approximate 
GCD computation with quadratic complexity w.r.t. polynomial degrees. 

1 Introduction 

The approximate polynomial greatest common divisor (denoted as AGCD) is 
a central object of symbolic-numeric computation. The main difficulty of the 
problem comes from the fact that is no universal notion of AGCD. One can 
find different approaches and different notions for AGCD. We will not give a 
review of all the existing work on this subject, but we will recall one of the most 
popular approaches to show how our work brings a different point of view on 
the problem. 

The main approach to the computation of an AGCD consists in considering 
two univariate polynomials whose coefficients are known with uncertainty. This 
uncertainty can be the result of the fact that the polynomials have floating point 
coefficients coming from previous computation (and so are subject to round-off 
errors). The most frequently adopted formulation is related to semi- algebraic 
optimization : given / and g two approximate polynomials, find two polynomials 
/ and g such that ||/ — /|| and — 5|| are small (lower than a given tolerance for 
instance) and such that the degree of gcd(/, g) is maximal. That is, one looks 
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for the most singular system close to the mput {f,g)- An e-gcd is obtained 
if the conditions ||/ — /|| < e and ||5 — 5|| < e are satisfied. One can try to 
compute the tolerance on the perturbation of the input polynomial thanks to 
direct computation (for instance from a jump on singular values of particular 
matrix for instance). This last approach has received a great interest following 
the work of Zeng using Sylvester like matrices f[14|). 

Here, we consider a slightly different problem. One of the polynomials, say 
/, is known exactly (it is the result of an exact model) and the second one, say g, 
is an approximate polynomial (result of measures or previous approximation for 
instance). This case occurs in applications such as model checking (to compare 
results of an exact model and measures). There are many other instances of 
such a problem, such as simplification of fractions when one of the polynomial 
is known exactly but the other one is not. 

We give a example of such a situation. When modeling an electromagnetic 
filter, one might want to parametrize its behavior with respect to the frequency. 
But one may need to do so even if there are singularities and to do so one may 
use Pade approximations of the electromagnetic signal at each point as a func- 
tion of the frequency. In some cases of interest, one can know all the singularities 
and so compute an exact polynomial called characteristic. Pade approximations 
are computed independently for each point by a numerical process and denom- 
inators may have a non trivial gcd with the "characteristic" polynomial. The 
denominators are not known exactly. So, in order to identify unwanted common 
factors in denominators one has to compute approximate gcds between an exact 
and non exact polynomials. 

This AGCD problem can also be interpreted as an optimization problem. 
Given / exactly and g approximately, compute a polynomial g close to g such 
that g has a maximal degree gcd with /. Our approach takes advantage of 
the asymmetry of the problem and of the structure of the quotient algebra 
C[x]/{f{x)) (more accurately, of the displacement rank of the multiplication 
operator in this algebra). So, we address the following problem : 

Problem 1 Let f{x) G C[x] a given polynomial and g{x) another polynomial. Find 
g{x) close to g{x) (in a sense that will be explained) such that f{x) and g{x) have a 
gcd of maximal degree. 

This may be also an interesting approach when one has two polynomials, one 
known with high confidence and another with worse accuracy. This approach 
may take advantage of this asymmetry which would not be possible for classical 
framework based on Sylvester or Bezout matrices. 

In this paper, we propose an approach and an algorithm to address this 
problem. The proposed algorithm is "fast" since the exponent of its complexity 
is better than the classical linear algebra exponent in the degree of the input 
polynomials. 

Organisation of the paper: The second section is devoted to some basic re- 
sult on algebra needed after, the third section gives an algebraic method for gcd 
based on linear algebra, the fourth section recalls the Barnett formula allowing 
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to compute the multiplication matrix without division, the fifth gives the dis- 
placement rank structure of the multiplication matrix, the sixth describes the 
final algorithm and experiments before finishing with conclusions and perpec- 
tives. 

2 Euclidian structure and quotient algebra 

In this section, we recall basic algebraic results allowing to understand the 
principle of our approach. All material in this section can be found (even in the 
non reduced case and in the multivariate setting) in 

Assume that K is an algebraically closed field (here we think about C). 

d 

Let f{x) and g{x) £ IK[x] and assume that f{x) ~ fd* 11 ^ d) ^^^d that 

Q ^ Q for aU i ^ j in {!,..., d}. Let A = K[x]/{f ) and tt : K[x] — > A be 

0(^-0 ) 

the natural projection. For i E {1, . . . ,d}, we define Li{x) ~ n'(C -C ) ' 

Lagrange polynomial associated to {^i, . ■ . , Cd}- Clearly, since deg(Li) < deg(/) 
we have 7r(Li) = L,, for all i G {1, . • . Let A* = HomK(A, K) be the usual 
dual space of A. For alH e {!,..., d}, we define l,;. : A — > K by lci(p) = 
p(Ci) for all p e A. The following lemma is obvious form the definition of the 

polynomials Li that for i and j G {1, . . . we have Li{Q) — 

This implies that the set {Li, . . . , L^} is a basis of A. A well known fact is 
that the set {1^^ , . . . , 1^^} form a basis A* dual of the basis {Li, . . . , Ld} of A. 
As a corollary, we have the Lagrange interpolation formula : Each p E A can 

d 

be written p{x) — X^IgCp) * Li{x). A funny consequence is that if we choose 

i=l 

{Li, . . . , Ld} as a basis of A, for all g E K[x], the remainder TT{g) of the euclidian 
division of 5 by / is given by (^(Ci), • ■ • , g{C,d)) in the basis {Li, . . . , Ld}, i.e. 

d 

r = 'Y^g(C,i)Li{x). In other word, divide (7 by / is equivalent to evaluate g at 

i=l 

the roots of /. 

The general philosophy of this last proposition will allows us to make a lot 
of proof in a very simple way. For example, it is very easy to see the different 
operation in A using this representation. Let g and h be to elements in A, then 

d d 

wehaveg + /i = Y^igid) + KCi))* Li{x) a.iidg*h = J2 (did) *HCi))Li{x) in A. 

4=1 i=0 

This allows us to avoid the use of the section a. In fact, the Lagrange polyno- 
mials Li, . . . ,Ld reveal a deeper structure on the algebra A : The polynomials 

Li, . . . , Ld are the idempotents of A, i.e. Li * Lj = 

Thanks to this description of the quotient algebra, it is easy to derive al- 
gorithms for both polynomial solving and gcd computation even though the 
problems are of very different nature. 



1 if i = j 
Oelse 



r Liiii = j 
\ Oelse 
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Remark that we have expressed everything in the monomial basis since it is 
the most widely used basis to express polynomials but we could use other bases. 
A particular basis is the Chebyshev basis where all results are exactly the same 
since it is a graduated basis. 

3 An algebraic algorithm for gcd computation 

To first give an idea on how to exploit the section above in order to design 
algorithm for gcd, We recall a classical method for polynomial solving (see |3] 
for instance). Proofs are given for the sake of completeness and because very 
similar ideas will lead us to the AGCD computation. 

3.1 Roots via eigenvalues 

d 

Let f{x) = J2fi^^ ^ ^ polynomial of degree d. Then we consider the 

i=Q 

matrix of the multiplication by x in C[x]/{f). Its matrix in the monomial basis 
1, . . . , x'^~^ is the following: 



Frob(/) 



( ^ ] 
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fd 



fd-1 I 

well known as the Froebenius companion matrix associated to /. 

Proposition 1 Let j[x) £ C[x]&e •polynomial of degree d with d distinct roots 
Z[f) — {zi,...,Zd}, then the eigenvalues o/Frob(/) are the roots of f{x), i.e. 
Spec(Frob(/)) = {zi, . . . , Zd} . 

Proof It follows directly from the fact that Frob(/) is the matrix of the mul- 
tiplication by x in C[x]/{f). But here we propose to give a direct proof by 
induction. In fact, we prove by induction that the characteristic polynomial of 
Frob(/) is f{x) itself (up to a sign and a scalar factor i.e.: 



1 

-X 1 





fd 
fd 



fd~l 
fd 
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since we have: 



-lo 
fd 

A 

/d 



/d-1 

fd 



/(x)-/o 



-X 1 







fd 



fd-1 
fd 



Jo 

fd 



and finally we have the wanted result. Re- 



-X 1 ••• 
-X 1 ■■■ 

••• -X 

and by assumption f{x) 
mark that this proof allows to avoid the condition that all the roots of f{x) are 
distinct. □ 

Then, to compute the roots of f{x) one can compute the eigenvalues of 
is Froebenius companion matrix. This is the object of the method proposed 
(reintroduced) by Edelman and Murakami [5] and revisited by Fortune [5] and 
many others trying to use the displacement structure of the companion matrix. 
In fact, often, the author realized that the monomial basis of the quotient algebra 
is not the most suitable one and proposed to express the matrix of the same 
linear application but in other basis. In the case of the Chebyshev basis this 
algorithm was already known by Barnett [2] and Cardinal later [3] . 

In the next section, we will also take advantage of the structure of the quo- 
tient algebra to design an algorithm for gcd computation mainly using linear 
algebra (eigenvalues are used in theory and never computed). 



-(a;*/(a;) + ^) 

Jd 



3.2 Structure of quotient and gcd 

Let f{x) andg{x) G K[x] such that they are both monic. As above, we denote 
A = K[a;]/(/) and d = deg(/). We denote denote {Ci, . ■ . , Q} the set of roots 
of f{x) and we assume that f{x) is squarefree, i.e. Q ^ Cj if * ^ j- We define 



Mg : I ^ ^ ^ TT{gh) ^^^^'^ "^^P) ^ ^ denote the remainder of p{x) G K[a;] 

by division by f{x). We denote Mg the matrix oi M.g in the monomial basis 
1, cc, . . . , a:''"^ of A but other bases can be used. A matrix representing the map 
Aig ]s called an extended companion matrix. 

Proposition 2 The eigenvalues of M.g are {.9(Ci): • • • j.9(Cd)}- 

Proof It is a direct corollary of the proposition ?? since if we write the matrix 
of this linear map in the Lagrange basis associated to {Ci, . . . , C^} is 

/ ff(Ci) ••■ 



V ■•■ 5(Cd) 
and gives the wanted result. □ 
Trivially, we have: 
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Corollary 1 We have corank(A^g) = deg(/) — rank(Alg) = deg(gcd(/, 5)). 

The column of index i of Mg is the column vector of the coefficients of 

x^~^ * g{x). 

Let pi, . . . ,pi be a basis of Ker(Mg) and let Pi{x),. . . , Pi{x) be the corre- 
sponding polynomials. First remark that AnnA(g') = {P{x) G K\P{x)*g{x) = 0} 
is an ideal of A. 

Lemma 1 The ideal AnnA(g') is a principal ideal. 
Proof Let us define 

s{x)= n (x-c). 

Cez(f)\(z(f)nz(g)) 

For all h e AnnA(fif) it is clear that Z{h) D Z{f)\{Z{f) n Z{g)) and then s 
divide h. Furthermore s £ AnnA(g') since in the Lagrange basis 

d 

s{x) * g{x) = ^s(Ci) * g{Ci)Li{x) = 0. 

i=l 

This shows that AnnA(g) — {s). □ 

To compute s{x), we built the matrix with columns formed hy pi, ... ,pi and 
we make a triangulation operating only on the columns. This way we obtain 
the polynomial of minimal degree linear combination of Pi (x), ... , Pi{x) and it 
is easily seen that this s{x) up to a multiplicative scalar factor. 

Lemma 2 The first column of a column echelon form of the matrix Kg built 
from a basis ofKei{Mg) is the generator of AuuA{g), i.e. it is the vector of the 
coefficients of s{x) up to a scalar multiplication. 

Proof Since the columns of a column echelon form of the matrix Kg are 
linearly independent, they form a basis of AnnA(g) as K-vcctor space. So s(x) 
is a linear combination of the polynomials associated to those columns. The 
polynomial associated to the column echelon form of Kg have all difFcrent degree 
(because it is an echelon form) and so s{x) is a linear combination of those 
polynomial. Because s{x) as the lowest degree possible, it is a scalar multiple 
of the polynomial associated to the first column. □ 

Proposition 3 f{x)Ag{x) = 

Proof By construction, we have s{x) * g{x) = Omod/(x) and so s{x) divide 
f{x). We also have gcd(|^, .g(.T)) = gcd(/(.x), ,g(x)) since the roots of are 

the root of f{x) where g{x) vanishes. Since deg(^|^) = deg{gcd{f{x), g{x)) we 
have the wanted result. □ 
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In all this section, we did not care if the polynomials are known in monomial 
or Chebyshev basis for instance. In fact, in order to have an algebraic algorithm, 
we only need to be able to perform euclidian division and this is always the case 
if the polynomial basis is graduated (as for monomial, Chebyshev, most of the 
orthogonal bases) . 

4 Bezoutian and Barnett's formula 

A classical matricial formulation of resultant is given by the Bezout matrix. 
In this part, we recall the construction of the Bezout matrix and a special 
factorization of the multiplication matrix expressed in the monomial basis. This 
factorization is called Barnett formula (see [2]). The Barnett's formula allows to 
build the classical extended companion matrix without using euclidian division 
and only stable numerical computations. Furthermore, this factorization reveals 
that the extended companion matrix has a special rank structure and we will 
use this fact later to design a fast algorithm to compute AGCD. 

Definition 1 Let f and g G C[a;] of degree m and n respectively (with n ^ m), 
we denote efjx,y) = /(^)gfa):/(^)g(^) ^ ^^.^.x^y^' = "e ^f,gA^)y' ■ The 
Bezout matrix associated with f and g is Bf^g ~ ( 9jj )ijg{o m-i}' 

Remark that since O f.g{x, y) = Q f,g{y, x) the matrix -B/,g is symmetric. The 
polynomials i^f,g,j(x) are univariate polynomials of degree at most m—1. One 
particular case of interest is when / = 1. In this case the Bezout matrix has a 
Hankel structure, i.e. 6*^^ — Oi^i^j^i. In this case we denote Hg^i{x) — Ki^g^i(x) 
for i £ {0, . . . , m — 1} which are called the Horner polynomials. 

Proposition 4 Let i G {0, . . . , m — 1}, the polynomial Hg^i{x) — ci,„i-i + • ■ • + 
Ci,ma;* has degree i and since they have different degree, they form a basis of 

m—1 

C[x]/{g). Furthermore, 1 ^g{x,y) ^ J2 Hg,m-i{x)y\ 

1=0 

Corollary 2 The matrix Bi g is the basis conversion from the Horner basis 
Hq, . . . , Mm~i to the monomial basis 1, x, . . . , x"^"'^ ofC[x]/{g). 

This leads us to the following theorem, known as Barnett formula (see [5]): 

Theorem 1 Let Mf be the multiplication matrix associated to f in C[x]/{g) in 
the monomial basis, we have: 

Mf = Bf^gBll. 

Proof Wehavee/,g(x,2/) = f{x)Siyl-^+g{x)^^^^^^ and so /(x)^^^!^ = 
Qf^g{x,y) inC[x,y]/ {g{x)). So, for each i G {0, . . . , m — 1}, we have &f^g^i{x) = 
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f{x)Qi^g^i{x). This last equality means that Bf^g is the matrix of the multipli- 
cation by f{x) in C[x]/{g). The result follows directly from this fact. □ 

The Barnett's formula reveals the rank structure of the multiplication ma- 
trix. Furthermore, this formula is already known if we choose Chebyshev basis 
instead of monomial basis to express the polynomials and the matrices have 
exactly the same nature. 

5 Structured matrices and asymptotically fast 
algorithms 

In this section, we briefly recall some basics on displacement structured matrices 
and related algorithms. 

5.1 Displacement structure 

Given an integer n and a complex number i9 with |t?| ~ 1, define the circulant 
matrix 

/ ^ \ 

1 



1 



V 10/ 

Next, define the Toeplitz-like displacement operator as the hnear operator 

. ^mxn ^ ^mxn 

Vt(A) = Zl^A - AZ^. 

A matrix A £ C™^" is said to be Toeplitz-like if Vt{A) is a small rank matrix 
(where "small" means small with respect to the matrix size). The number 
a = rank(V(A)) is called the displacement rank oi A. If A is Toeplitz-like, then 
there exist (non- unique) displacement generators G G c™x" q^^lA H £ c^xn 
such that 

V{A) = GH. 

Toeplitz matrices and their inverses are examples of Toeplitz-like matrices. An- 
other useful example is the multiplication matrix Mf, which has Toeplitz-like 
displacement rank equal to 2, regardless of its size. 

A similar definition holds for Cauchy-like structure; here the relevant dis- 
placement operator is 

Vc(A) = DiA~AD2, 

where Di and D2 are diagonal matrices of appropriate size with disjoint spectra. 
See [lO] for a detailed description of displacement structure. 
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5.2 Fast solution of displacement structured linear sys- 
tems 

Gaussian elimination with partial pivoting (GEPP) is a well-known and reli- 
able algorithm that computes the solution of a linear system. Its arithmetic 
complexity for an n x n matrix is asymptotically 0{n^). But if the system ma- 
trix exhibits displacement structure, it is possible to apply a variant of GEPP 
with complexity 0{n^). The main idea consists in operating on displacement 
generators rather than on the whole matrix; see [7] for details. 

Strictly speaking, the GKO algorithm performs GEPP (or, equivalently, 
computes the PLU factorization) for Cauchy-like matrices. However, several 
authors have pointed out (see [7], [S], [H]) that Toeplitz-like matrices can be 
stably and cheaply transformed into Cauchy-like matrices; the same is true for 
displacement generators. 

Consider, for instance, the case = 1 and let ^ be an n x n Toeplitz-like ma- 
trix with generators G and H. Denote by Dq the matrix diag(l, e^'/", . . . , e("~i)'^*/") 
and let F be the Fourier matrix of size n x n. Then the matrix FAD^^F^is 
Cauchy-like, of the same displacement rank as A, with respect to the dis- 
placement operator defined by Di = diag(l, e^'^'/"-, . . . , e^'^^C"-!)/") and D2 = 
diag(e''*/", e^'"'/", . . . , e^^""^)'"*/"). Its Cauchy-like generators can be computed 

as G = FG and h" = FDqH" . 

Generalization to the case oi m x n rectangular matrices is possible. In this 
case, the parameter -d should be chosen so that the spectra of Di and D2 are 
well separated (see [1] and [3]). 

We also point out that the GKO algorithm can be adapted to pivoting tech- 
niques other than partial pivoting ([8], [E]). This is especially useful in case of 
instability due to internal growth of generator entries. A Matlab implementa- 
tion of the GKO algorithm that takes into account several pivoting strategies is 
found in the package DRSolve described in [I]. In our implementation, we use 
the pivoting strategy proposed in [5]. 

6 A structured approach to AGCD computation 

We propose here an algorithm that exploits the algebraic and displacement 
structure of the multiplication matrix to compute the AGCD of two given poly- 
nomials with real coefficients (as defined in section [T|). 

6.1 Rank estimation 

It has been pointed out in Section [3] that the rank deficiency of the multiplica- 
tion matrix equals the AGCD degree. Here we use the structured pivoted LU 
decomposition to estimate the approximate rank of the multiplication matrix. 
Recall that Mg has a Toeplitz-like structure with displacement rank 2; it can 
then be transformed into a Cauchy-like matrix Mg as described in Section 5.2. 
Fast pivoted Gauss elimination yields a factorization Mg = P1LUP2, where L 
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is a square, nonsingular, lower triangular matrix with diagonal entries equal to 
1, U is upper triangular and Pi, 2 are permutation matrices. Inspection of the 
diagonal entries (or of the row norms) of U allows to estimate the approximate 
rank of Mg and, therefore, of Mg. 



6.2 Minimization of a quadratic functional 

Let us suppose that: 

• the polynomial f{x) — X]j=o fj^^ exactly known, 

• the polynomial g{x) = jyjLogjX^ is approximately known and may be 
perturbed, so that we consider its coefficients as variables, 

• the AGCD degree is known. 

Then we can reformulate the problem of AGCD computation as the minimiza- 
tion of a quadratic functional. Indeed, recall that the cofactor v{x) with respect 
to f{x) is defined by the "shortest" vector (i.e., the vector with the maximum 
number of trailing zeros) that belongs to the null space of Mg. We assume v{x) 
to be monic; we denote its degree as k and we have 



MgV = Mg 



( 



( \ 



Vk-l 

1 




V J 



V / 



Also observe that the entries of Mg are linear functions of the coefficients of 
g{x). Then the equation MgV = can be rewritten as J^{g,v)=0, where the 
functional is defined as 

T : C'"+^ X C'' — >R+ 
T{g,v)^\\Mgv\\l 

For a preliminary study of the problem, we have chosen to solve the equation 
J^{g, v)—0 by means of Newton's method, applied so as to exploit structure. 
Denote by z = [go, ■ ■ ■ , gm,Vo, ■ ■ ■ ,Vk-i]'^ the vector of unknowns; then each 
Newton step has the form 



Jig 



0) 



In particular, notice that the Jacobian matrix associated with is an n x (m + 
A; + 1) Toeplitz-like matrix of displacement rank 3. This property allows to 
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compute a solution of the linear system J{g'^^\v^^^)y = Mg(j)V^J^ in a fast way; 
therefore, the arithmetic complexity of each iteration is quadratic w.r.t. the 
degree of the input polynomials. 

We propose in the future to take into consideration other optimization meth- 
ods in the quasi-Newton family, such as BFGS. 

6.3 Computation of displacement generators 

In order to perform fast factorization of the multiplication matrix Mg, we need 

to compute Toeplitz-like displacement generators. It turns out that the range 
of V{Mg) is spanned by the first and last column of the displaced matrix, and 
the columns of indices from 2 to n — 1 are multiples of the first one. Therefore, 
it suffices to compute a few rows and columns of Mg in order to obtain displace- 
ment generators. This can be done in a fast and stable way by using Barnett's 
formula. If we denote as Cj the jth vector of the canonical basis of C , then 
the computation of the 7-th column of M^can be seen as 

Mg(:,j)=S(/,,9)-(B(l,/)-ie,), 
that is, it consists in solving a triangular Hankel linear system and computing 
a matrix- vector product. For row computation, recall that the Bezoutian is a 
symmetric matrix; we have analogously: 

M,(i,:) =eJ-B(/,5)-B(l,/)-^ = {B{lJ)-^B{f,g)e'^f , 
so that the computation of a row of Mg amounts to performing a matrix-vector 
product and solving a Hankel triangular system. 

A similar approach holds for computation of displacement generators of the 
Jacobian matrix J{g,v) associated with the functional J^{g,v). 

6.4 Description of the algorithm 

Input: coefficients of polynomials f{x) and g{x). 
Output: a perturbed polynomial g{x) such that / and g have a nontrivial com- 
mon factor. 

1. Estimate the approximate rank k of Mg by computing a fast pivoted LU 
decomposition of the associated Cauchy-like matrix. 

2. Again by using fast LU, compute a vector v = [vo,vi, , Wfe_i, 1, 0, , 0]^ 

in the approximate null space of Mg. 

3. Apply structured Newton with initial guess {g, v) and compute polynomials 
g and v such that / and g have a common factor of degree deg f — k and v 
is the monic cofactor for /. 

6.5 Numerical experiments and computational issues 

We have written a preliminary implementation of the; proposed method in Mat- 
lab (available at the URL http : //www.unilim. f r/pages4)erso/paola.boito/MMgcd.m). 
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The results of a few numerical experiments are shown below. The polynomi- 
als / and g are monic and have random coefficients uniformly distributed over 
[—1,1]. They have an exact GCD of prescribed degree. A perturbation is then 
added to g. The perturbation vector has random entries uniformly distributed 
over [— ?7, rj] and its norm is of the order of magnitude of rj. We show: 

• the residual J'{g,v), 

• the 2-norm distance between the exact and the computed cofactor v, 

• the 2-norm distance bcrtwecm the exact and the computed perturbed poly- 
nomial g (which is expected to be roughly of the same order of magnitude 
as T]). 

In the following table we have taken r] =le-5. 



n, TO, deg / 


^{9, v) 


\\v-i\\2 


11.9 -.9II2 


8, 7,3 


1.02e-15 


1.19e-15 


1.40e-5 


15, 14, 5 


1.51e-15 


2.26e-15 


1.35e-4 


22, 22, 7 


2.07C-13 


1.40e-13 


2.20e-4 


36, 36, 11 


1.19e-12 


5.07e-14 


0.0012 


Here are results for rj =le-8: 


n, TO, deg / 




\\v-i\\2 


11.9 -.9II2 


8, 7,3 


5.49e-15 


1.63e-15 


5.85e-8 


28, 27, 13 


7.90e-14 


8.98e-14 


6.50e-7 


38, 37, 13 


4.88C-12 


4.26C-12 


2.30C-5 


58. 57, 23 


2.03C-12 


4. 40c- 12 


2.54C-4 



There are several issues in our approach that deserve further investigation. 
Let us mention in particular: 

• The choice of a threshold (or a more refined technique) for estimating 
approximate rank. 

• Normalization of polynomials: here wc mostly work with monic polyno- 
mials, but other normalizations may be considered. 

• The structured implementation of the optimization step (minimizing J^{g, v)). 
We have used for now a heuristic structured version of the Gauss- Newton 
algorithm. Observe that each step of classical Gauss-Newton applied to 
our problem has the form z^^'^^'^ = z^^^ —y'^^\ where z^^^ is the vector con- 
taining the coefficients of the j-th iterate polynomials g'^^^ and f ^■'^ and y^^^ 

is the least-norm solution to the underdetermined system J{g^-''^ , f ^■'^ )y^^^ = 
MgU) v^^^ . Computing this least-norm solution in a structured and fast way 
is a difficult point that will require more work. Our implementation gives 
a solution which is not, in general, the least-norm one, even though it 
is typically quite close. Further work will also include a study of other 
possible optimization methods that lend themselves well to a structured 
approach. 
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7 Conclusions 



We have proposed and implemented a fast structured matrix-based approach to 
a variant of the AGCD problem, namely, the problem of computing an approx- 
imate greatest common divisor of two univariate polynomials, one of which is 
known to be exact. To our knowledge, this variant has been so far neglected in 
the existing literature. It may be also interesting when one polynomial is known 
with high accuracy and the other is not. 

Our approach is based on the structure of the multiplication matrix and 
on the subsequent reformulation of the problem as the minimization of a siiit- 
ably defined functional. Our choice of the multiplication matrix Mg over other 
resultant matrices (e.g., Sylvester, Bezout...) is motivated by 

• the smaller size of Mg, with respect e.g. to the Sylvester matrix, 

• the strong link between the null space of Mg and the gcd, and in particular 
the fact that the null space of Mg immediately yields a gcd cofactor, 

• the displacement structure of Mg, 

• the possibility of computing selected rows and columns of Mg in a stable 
and cheap way, thanks to Barnett's formula. 

This is, however, a preliminary study. Fiirther work will include generalizations 
of the proposed problem and a more thorough analysis of the optimization part 
of the algorithm. Furthermore, this approach can be generalized in several 
intersting way: 

• using better bases then the monomial one, 

• it can be extended to some multivariate setting to compute the co-factor 
of a polynomial g in C[xi, . . . , . . . , /„) when /i, . . . , /„ define a 

complete intersection since Barnett formula still holds, 

• to compute the AGCD of / with gi, ■ ■ ■ ,gk where / is known with accuracy 
but gi,- ■ .,gk are inaccurate, one can take 5 as a linear combination of 
gi,...,gk with our method and succeed with a high probability. 
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